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We report the bulk and surface properties of lithium computed within a full potential LCGTO 
formalism using both density functional theory and the Hartree-Fock approximation. We examine 
the convergence of computed properties with respect to numerical approximations and also explore 
' the use of finite temperature density functional theory. We demonstrate that fully converged calcu- 

lations reproduce cohesive properties, elastic constants, band structure, and surface energies in full 
agreement with experimental data and, where available, previous calculations. 



<3\ 



3 



O 

CZ3 



-a 



I. INTRODUCTION 



Lithium has been the subject of considerable interest 
over many years. Although its electronic structure is rel- 
atively simple, its structural properties still pose a sienifi,- 
cant challenge to both experimentuTa and simulationDtJ. 
Lithium is very soft; the determination of its elastic con- 
stants and surface energies requires experiments of high 
accuracy and simulations of high numerical stability. The 
calculation of surface formation energies is particularly 
delicate as has been discussed recentlyE3E3. 
q ■ The aim of this article is to present a comprehensive 
and systematic study of the band structure, cohesive en- 
ergy, elastic constants, phase stability, and surface ener- 
gies of lithium. Very few systematic studies of the de- 
pendence of results on the computational parameters are 
available. This however is especially important when en- 
ergy differences are required as for example in calculation 
r*"** ' of surface energies or phase stabilities. We present the 
results of fully converged, full potential, all-electron cal- 
culations based on both density functional theory and the 
Hartree-Fock approximation. We examine the use of fi- 
nite temperature density functional theory as a technique 
for accelerating convergence with respect to reciprocal 
""j" 1 , space sampling. We expand the crystalline orbitals as a 
linear combination of Gaussian type orbitals (LCGXOl 
This approach is very well established for insulatorsc£l — . 
We find, in accord with a recent study of magnesiumEJ, 
that this approach is also well suited to the simulation of 
a free-electron metal. 

The paper is organized as follows: in section 0, we 
discuss the computational parameters. In sections |Tn 



and |iyj, we discuss results for lithium bulk and surfaces, 
respectively, and summarize our conclusions in section [v|. 



II. BASIS SET AND METHOD 



All the calculations were performed with the program 
package CRYSTAiO. The main numerical approximation 
in our approach is the choice of the Gaussian basis set. 
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The difficulties of selecting basis sets for .metallic systems 
have been explored in previous studiesE3E3. In princi- 
ple the quality of a calculation can be systematically im- 
proved by adding additional functions to the basis set and 
optimizing their exponents in a suitable reference state - 
usually the bulk crystal. In practice one must balance the 
overcomplctcness of the basis set, which leads to linear 
dependence, with the need for additional variation free- 
dom. For molecular systems and insulating solids these 
problems have largely been overcome and Gaussian basis 
sets are very widely used. For metallic systems and their 
surfaces in particular there have been very few systematic 
studies. 

In a solid the tails of the atom centered functions over- 
lap strongly and so diffuse basis functions optimized for 
the description of atomic or molecular systems are not 
useful and indeed may give rise to linear dependence. 
We are thus unable to simply use basis sets from the 
many libraries developed for the description of molecular 
systems. We have therefore developed an hierarchy of 
basis sets of increasing quality in order to examine the 
convergence of computed properties. 

The smallest basis set used has 3 s-symmetry functions 
and 2 p-symmetry functions and is denoted as [3s2p| . 
The [Is] radial function was taken from Ref. E4 The 
exponents of the two additional sp-shells were optimized 
in local density appropdmation (LDA) calculations (with 
Dirac-Slater exchanges and Perdew-Zunger's correlation 
functional^) for the solid at the experimental lattice con- 
stant. The lowest energy was obtained with exponents 
of 0.50 and 0.08. Similar results were obtained when us- 
ing the Per-dew Wang gradient corrected approximation 
(PWGGA)eI. However, as an exponent of 0.08 gives rise 
to a very diffuse basis function close to numerical insta- 
bility, instead exponents of 0.50 and 0.10 were chosen. 
This [3s2p] basis set is very robust and computationally 
efficient — it does not give rise to linear dependence even 
when the bulk is strongly distorted (for example, to de- 
termine the elastic constants). 

A [4s3p] basis set was obtained by using three expo- 
nents (0.50, 0.20, and 0.08) - which were chosen to be 
"even tempered", i.e. the ratio between the exponents 
is kept fixed (2.5 in this case). This ratio is close to 
the lowest which can be tolerated before on-site (atomic) 
linear dependence is seen. It is however also known to 
converge the atomic energy to within less than 10~ 4 Eh 
(Eh = 27.2114 eV) of the exact Hartree-Fock ground 
state energy (see the analysis in Ref. |27]) . Finally, an ad- 
ditional polarization function of d-symmetry was added 
and the exponent optimized within a PWGGA calcula- 
tion to be 0.15. However, the d-function leads only to a 
minor change in the total energy. The energy varies only 
by 5 x 10~ 5 i?^ when, e.g., changing the exponent to 0.5. 
The basis sets developed in this manner are displayed in 
table g. _ 

Both at the Hartree-Fock (HF) and B3LYPEI (invov- 
ing a hybrid of Fock exchange and a modificatieia-pf the 
Becke gradient corrected exchange functiona]ESE2l, and 
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the Vosko-Wilk-Nusair local correlation functional veil 
and the gradiept corrected correlation potential by Lee, 
Yang and ParrE3) levels, an optimization of basis set ex- 
ponents was not possible. Instead, the outermost ex- 
ponent became more and more diffuse until finally the 
solution became unstable. This is a well known pathol- 
ogy of the use of Fock exchange in metallic systems (see 
also the discussion in Ref . ^3|) . When features of the HF 
solution are discussed in this article, they were obtained 
with the [3s2p] basis set (outermost exponents 0.50 and 
0.10). 

In order to compute binding energies the free atom 
is calculated within a spin-polarized formalism with the 
same [Is] function but with additional s-exponents 0.60, 
0.24, 0.096, 0.04, and 0.016 to describe the long range 
behaviour of the atomic wavefunction. 

For the LDA and PWGGA calculations we also expand 
the exchange and correlation potentials in an auxiliary 
Gaussian basis set which consists of 13 even tempered 
s-functions with exponents from 0.1 to 2000, 3 even tem- 
pered p-functions with exponents from 0.1 to 0.8, and 2 
c?-functions with exponents of 0.12 and 0.3. This is suf- 
ficient to integrate the charge density to an accuracy of 
10 -7 |e|. For the free atom, we use an auxiliary basis set 
with 18 even tempered s-functions with exponents from 
0.0037 to 4565. 

Reciprocal space sampling is a delicate problem espe- 
cially in metals— - The sampling is performed on a Pack- 
Monkhorst netcJ where the density of points is deter- 
mined by a shrinking factor. The Fermi energy and 
shape of the Fermi surface are determined by interpola- 
tion onto a "Gilat" net. This net is simply related to the 
Pack-Monkhorst net by an additional subdivision factor. 
To further improve convergence, the finite-,temperature 
generalization of density functional theoryEJ can be used 
to apply Fermi surface smearingpj. In table Q, the de- 
pendence of the total energy on the density of points in 
the Pack-Monkhorst net is displayed. For this purpose, 
PWGGA calculations on a body centered cubic (bcc) lat- 
tice at the equilibrium lattice constant of 3.44 A were 
performed. 

At zero temperature, even with the largest net used, 
the energy is still slightly decreasing as more points are 
used. For a smearing of 0.001 Eh (which corresponds to 
a temperature of T = 0,0 "* g " = 316K with Boltzmann's 

constant fcs = 3.1667 x 10~ 6 -^), the energy is stable up 
to a few /iEh- A higher number of sampling points in 
the Gilat net leads to a systematic improvement at zero 
temperature. At finite temperature, the number of sam- 
pling points in the Gilat net does not influence the results 
when a sufficiently high number in the Pack-Monkhorst 
net is chosen. As shown in table 0, the difference in en- 
ergy for the different number of sampling points in the 
Gilat net is of the order of only a few \iEh for a fixed 
shrinking factor of 24 in the Pack-Monkhorst net. Note 
that for the most dense net at UbT = 0.001-E/,,, the differ- 
ence in energy between smeared and unsmeared results 
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is less that 10~ 4 Eh- At an even higher temperature of 
0.02 Eh, the energy converges to a value which is 6 mEh 
higher than at 0.001 Eh- An estimate of the zero tem- 
perature energy is possible by using the approximation 
E(0K) = \(E(T) + F(T)) (with F = E - TS being 
the free energy and exploiting the fact that the energy 
increases quadratically with temperature for low temper- 
ature) as suggested in Ref. |(| The electronic entropy S 
is defined as 

S = k B fi In fi + (1 -fi) ln(l - fi) 

with fi being the Fermi function. This leads to a value 
of E(0) extrapolated from ksT = 0.02Eh which deviates 
by less than lO -3 ^ from the value obtained at a tem- 
perature of 0.001 E h . The functions E(T), F(T) and 
h(E(T) + F(T)) for a fixed value of 145 sampling points 
(corresponding to a shrinking factor of 16 in the Pack- 
Monkhorst net) are also displayed in Figure [|. Indeed, 
even for relatively high temperature, E(0) is well approx- 
imated by \{E{T) + F(T)). 




' 0.000 0.020 0.040 0.060 

Temperature [E h ] 

FIG. 1. Energy E(T), free en- 

ergy F(T) and \{E{T) + F{T)) of Li bulk, bcc lattice. A 
shrinking factor of 16 was used. 

In conclusion, as we need a high accuracy for the en- 
ergy difference between the different phases of Li, we 
chose for the calculations on lithium bulk a shrinking 
factor of 16 for the Pack-Monkhorst net which gives 145 
sampling points in the irreducible Brillouin zone of the 
non-distorted bulk and a temperature of 0.001 Eh- This 
ensures that convergence of the energy to at least lO -4 ^ 
with respect to reciprocal space sampling is achieved. 



III. RESULTS FOR BULK LITHIUM 
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A. Band structure 



Li: LDA 
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FIG. 2. LDA band structure at the equilibrium lattice con- 
stant, bcc lattice 

In figure || the LDA band structure for the bcc 
structure is displayed. The occupied bands and the 
lower unoccupied bands are in excellent agreement with 
results earlies-, obtained with Gaussian basis sets and 
X a exchangeO, by the Kohn-Korringa-Rostoker (KKR) 
method^ and Slater exchange, augmented plane wave 
(APW)iJ and modified APW (MAPW)B calculations 
both using iDA, or linear muffin tin orbital (LMTO) 
calculations^] with a combination of exchange using the 
Langreth and Mehl functional^ and LDA correlation. 
When using PWGGA instead of LDA, the band struc- 
ture does not exhibit major differences. The_experimen- 
tally known conduction band width (~ 4 eV)e3 is slightly 
lower than the result calculated here (4.6 eV). The slow 
decay of the density matrix in metals leads to difficulties 
when Fock exchange is involved: the summation of the 
exchange series in direct space is very long ranged and 
is truncated at a large but finite distance. This cutoff 
for large distances in direct space results in numerical in- 
stabilities for small k. Thus, both in the Hartree-Fock 
and B3LYP band structures, artificial oscillations can be 
found around the Gamma point. As usual for Hartree- 
Fock calculations, we find that the bandwidth is roughly 
twice as large as the experimental bandwidth. 



B. Cohesive properties 



Table III gives results for ground state properties of 
bcc Li (lattice constant, bulk modulus, cohesive energy 
and elastic constants Cn, C\\ —Cm, and C44). The elas- 
tic constants were obtained by applying a rhombohedral 
distortion for C44, a tetragonal distortion for Cn, and an 
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orthorhombic distortion for C\x — Cn to the solid at the 
equilibrium lattice constant. Cohesive energy and lat- 
tice constant agree well with experiment and depend only 
weakly on the basis set. The bulk modulus is obtained 
from the energy as a function of volume and agrees within 
the accuracy of the fit with that obtained using the rela- 
tion B = \{Cn + 2C12). Elastic constants have a strong 
dependence on the basis set and the deviation from ex- 
periment improves especially when going from [3s2p] to 
[4s3p]; the d- function has only minor impact on the re- 
sults. We find that PWGGA is closer to experiment with 
results similar to Ref. 0. 



C. Relative stabilities 

The experimental crystal structure of lithium at zero 
temperature is still unclear (at room temperature a bcc 
structure is favoured) . Both face-centered cubic (fee) and 
hexagonal (hex) structures have been suggested as well as 
more sophisticated structures such as 9R (a niriedayer se- 
quence of close-packed planes ABCBCACAB)c3 or mix- 
tures of these phases. Therefore, we also investigated the 
relative stability of the different phases. In Figure ||. to- 
tal energies of bcc, fee and hex phase, obtained with the 
PWGGA functional and the best basis set ([4s3pld]), are 
displayed. 



Li: PWGGA 

Relative Stability of the Phases 



■H- -7.5300 




x bcc 
+ fcc 
□ hex 



17.0 



19.0 21.0 
Volume/Atom [A 3 /atom] 



FIG. 3. Relative stabilities of the different phases as a func- 
tion of volume, [4s3pld] basis. 



The c/a ratio of the hexagonal phase remains close to 
the ideal close packed value of 1.633 (it varies between 
1.631 and 1.635 which is within the accuracy of the fit). 
We find that the closed-packed structures are slightly 
lower in energy than the bcc structure. Although it is 
reasonable to conclude that the close packed structures 
are lower in energy than bcc it is not possible to resolve 
the difference in energy between the fee and hex phases. 
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The variation of the energy with basis set is given in table 



IV . As the basis set is systematically improved the energy 
difference between bcc and fee increases from 1 xlO~ 5 Eh 
to 4 xlCP 5 Eh, and between bcc and hex decreases from 
7 xlCT 5 ^ to 4 xlCT 5 ^. We note that even the d- 
function influences this energy splitting. These results, 
which do not include the zero point energy, indicate a 
preference for closed-packed structures whichis-iii agree- 
ment with most of the previous calculationsotil except 
for oneEa (see table |v|). 

HF calculations were only possible with the smallest 
[3s2p] basis set where the order of phases is different with 
hep being the lowest in energy, followed by bcc and fee 
being highest. The same is found in LDA with the small- 
est [3s2p] basis set, but changes when the basis set is in- 
creased to [4s3p] where both fee and hep are 2 xlO~ 4 Eh 
lower in energy than bcc (again, calculations with the 
best [4s3pld] basis set were not possible because of lin- 



ear dependence) — table [IV 



As the energy splitting between the close packed and 
bcc phases is rather small the zero-point energy differ- 
ence cannot be neglected. The first published calculation 
of this, based on the harmonic approximation^ gave an 
additional stabilization of 9 x 10~ 5 Eh of the fee phase 
compared to the bcc phase. In these calculations how- 
ever the hep phase was found to be higher in energy that 



the bcc phase (table IV). Very recently, in a calculation 
also including anharmonic effects, the stabilization was 
calculated to be about 1.5 x 10 _5 £a, both for fee and 
9R phases relative to the bcc phasefl. In addition, the 
authors computed the variation of the vibrational free 
energy as a function of temperature and found a phase 
transition from a closed packed to a bcc phase at a tem- 
perature of T ~ 20QK. 



IV. RESULTS FOR SURFACES 

A further test of this approach is the calculation of 
surface energies. We model lithium surfaces by using a 
slab and varying the numbers of layers (a bcc structure 
is assumed). Surface energies can be calculated in two 
ways, either by deriving a bulk energy by subtracting 
energies of two slabs with n and to layers: 

1 Ti 

E S urface = -z{ E slabijl) - {E s i ab {n) - E sUlb {n - m)) — ) 

2 TO 

(1) 

which has the advantage of a systematic error cancella- 
tion (in particular the reciprocal space sampling is con- 
sistent between the bulk and slab energies) or by using 
an independent bulk energy 

Efface = l(E slab (n)-E bulk xn) (2) 

All the quantities E surface , E s i ab (n), and E bu i k are ener- 
gies per atom. 
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In Figure ]^ results for surface energies using both 
equation [l] (with m=l) and || are displayed. Equation 
|l| leads to relatively strong oscillations (dotted line with 
pluses) and gets more stable the larger to. Numerical 
noise in the expression (E s i ab (n) - E s i ab (n - m))~ is 
reduced for larger values of to. Equation however, 
leads at zero temperature to a slight linear decreasing of 
the surface energy as a function of the number of lay- 
ers (thin line without additional symbol). The reason 
for the nonvanishing slope is that the energy difference 
E(n) — E(n — 1) is not identical to the energy of the 
bulk due to the systematic errors in the convergence of 
the total energy with respect to reciprocal lattice sam- 
pling (this was also emphasized in Refs. |l6|,[l7]). As shown 
in table [n|, at zero temperature the bulk energy is still 
changing on the order of 10 Eh with increasing sam- 
pling point number. Similarly the bulk energy varies 
when extracted from the slab. This slight discrepancy 
gives rise to a variation of the surface energy with the 
number of layers with a slope of 10~ 4 Eh/ atom /layer. 

The origin of the poor convergence with respect to re- 
ciprocal lattice sampling is due to the sharp cutoff im- 
posed by the Fermi energy. One possibility to obtain 
the surface energy would be to use the intercept from 
a linear fit, but a better and simpler way to alleviate 
the difficulties associated with reciprocal lattice sam- 
pling is smearing the Fermi surface with a finite temper- 
ature. Already at a temperature of 0.001 Eh, the slope 
is virtually zero (thick line without additional symbol). 
This is consistent with table [H] as the bulk energy con- 
verges much faster at finite temperature. At a higher 
temperature of fc^T = 0.02Eh, E(T) clearly deviates 
from E(0) (thin line with stars) and the approximation 
E(0) = \{E{T) + F(T)) should be applied. This works 
very well when comparing the corresponding results (thin 
line with crosses) with results calculated at ksT = 0.001 
Eh (thick line without additional symbol). The surface 
energy obtained this way is only slightly higher than that 
from a calculation at fc^T = 0.001 Eh which is consistent 
with figure [l]. 

A higher smearing temperature also reduces the oscil- 
lations in the curves both when using equation [l] or |2[ It 
even allows to substantially reduce the number of sam- 
pling points as shown in Figure ^| At low temperature, 
a high number of sampling points is necessary to obtain 
the correct result whereas at high temperature already 
a shrinking factor of 4 (resulting in 6 sampling points in 
the irreducible zone) is sufficient. It should be noted that 
in this case we used Equation [l] so that all the data is 
consistently extracted from calculations on slabs. 

In conclusion, calculations at zero temperature are 
very cumbersome when results from calculations on slabs 
and on the bulk have to be combined. The error cancella- 
tion can be maximized by extracting the surface energy 
from calculations on slabs only. The problem of error 
cancellation between bulk and slab is already improved 
at very low temperature when it is possible to fully con- 
verge the bulk energy with respect to reciprocal lattice 
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sampling. Higher temperature also leads to a smoother 
behaviour of the surface energy as a function of number 
of layers. Finally, calculations at high temperature can 
be performed with a strong reduction of the number of 
sampling points and using E(0) = \{E(T) + F(T)) as an 
approximation for zero temperature results as suggested 
in Ref. g 

In table [v], LDA and PWGGA results for the unrelaxed 
(100), (110) and (111) surface are summarized. The lat- 
tice constant was chosen as the bulk equilibrium lattice 
constant, a temperature of 0.001 Eh and a higher shrink- 
ing factor of 24 resulting in 91 sampling points for the 
slabs and 413 sampling points for the bulk was used. 



Li (100) Surface 

PWGGA 



+ fusing E(n), E(n-1), k B T=0 E h 

using E(n), E(bulk). k B T=0 E„ 

using E(n), E(bulk), k B T=0.001 E„ 

x x using E(n), E(bulk), k B T=0.02 E h 

X X using (E(nHF(n))/g, (E(bulk)+F(bulk))/2, k n T=0.02 E„ 




10 20 30 

Number of Layers 

FIG. 4. (100) Lithium surface energy with a shrinking fac- 
tor of 24 in the Pack-Monkhorst net, [3s2p] basis. 



Li (100) Surface 

PWGGA 




0.02 
Temperature [EJ 



FIG. 5. (100) Lithium surface energy with two different 
reciprocal lattice samplings extracted from two slabs with 5 
and 6 layers using \{E(T) + F(T)) and Equation |l] 
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As found in Ref. ^ for the jellium model without 
additional long-range corrections, surface energies in 
PWGGA are lower than in LDA (the different lattice con- 
stants for PWGGA and LDA are not the reason, an eval- 
uation at the LDA equilibrium lattice constant leads to a 
change of the PWGGA surface energy which is negligible 
compared to the difference between LDA and PWGGA 
surface energies). The energies are reduced when going 
from the smaller [3s2p] to the [4s3p] basis set, introduc- 
ing a d function leads only to-sappr changes. Our results 
agree well with the literatureE3E3. 



V. CONCLUSION 

In this study of lithium metal, we presented accurate 
results using a full potential, all electron density func- 
tional scheme. Results for cohesive properties, elastic 
constants, band structure, and surface energies are in full 
agreement with experiment and calculated values from 
literature. The results in best agreement with exper- 
iment were obtained with the gradient corrected func- 
tional of Perdew and Wang. Hartree-Fock calculations 
for lithium are very difficult as it is impossible to opti- 
mize the exponents because of the very long range of the 
exchange interaction; the same problems appear in func- 
tional involving an admixture of Fock exchange. We 
demonstrated the convergence of the different properties 
with respect to the computational parameters by using 
an hierarchy of basis sets, different reciprocal lattice sam- 
plings and different smearing temperatures. We showed 
that finite temperature calculations can be used to im- 
prove convergence and still an extrapolation to zero tem- 
perature is possible and accurate. Quantities like cohe- 
sive energy and lattice constant are already stable with 
the smallest basis set; elastic constants and surface en- 
ergies are more sensitive. The most difficult quantity 
to calculate is the energy splitting between the differ- 
ent phases where we have reached the limit of numerical 
accuracy. We can not make a prediction about the pre- 
ferred crystal structure; the energy difference is so small 
that, from the computational point of view, even subtle 
changes such as introducing a d-function_.are important 
and zero point energies must be includecoQ. We confirm 
the finding of Ref. |2C] that an approach based on Gaus- 
sian type functions provides a reliable and very efficient 
description of metallic systems. 
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TABLE I. Basis sets 



exponent s contraction p contraction d contraction 

[Is] 840.0 0.00264 

217.5 0.00850 

72.3 0.0335 

19.66 0.1824 

5.044 0.6379 

1.5 1.0 

[3s2p] 

[2sp] 0.50 1.0 1.0 

[3sp] 0.10 1.0 1.0 

[4s3p] and [4s3pld] 

[2sp] 0.50 1.0 1.0 

[3sp] 0.20 1.0 1.0 

[Asp] 0.08 1.0 1.0 
[d] 0.15 1.0 
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TABLE II. Convergence of the total energy with respect to the number of sampling points. The results are obtained from 
PWGGA calculations with a [3s2p] basis set, for the bcc lattice at a lattice constant of 3.44 Aand for three temperatures Eh, 
0.001 E h , and 0.02 E h . 



shrinking number of Shrinking factor of Gilat net 

factor of the sampling points in in multiples of shrinking factor 

Pack Monkhorst irreducible part of the of the Pack-Monkhorst net 

net Pack-Monkhorst net 



k B T = 0E h ; E{T) = E{0) 







1 


2 


3 


4 


4 


8 


-7.520822 


-7.520213 


-7.520284 


-7.520373 


8 


29 


-7.527049 


-7.529923 


-7.530344 


-7.530488 


12 


72 


-7.528987 


-7.529867 


-7.530023 


-7.530075 


16 


145 


-7.529574 


-7.530032 


-7.530125 




18 


195 


-7.529697 


-7.530077 






20 


256 


-7.529796 


-7.530101 






24 


413 


-7.529915 


-7.530130 







k B T = 0.00lE h ; E{T) 



4 


8 


-7.524523 


-7.522697 


-7.521528 


-7.521029 


8 


29 


-7.529710 


-7.530632 


-7.530634 


-7.530654 


12 


72 


-7.529982 


-7.530156 


-7.530141 


-7.530128 


16 


145 


-7.530175 


-7.530183 


-7.530181 




18 


195 


-7.530195 


-7.530189 






20 


256 


-7.530185 


-7.530188 






24 


413 


-7.530200 


-7.530185 







k B T = 0.001£ h ; E(0) = ±(£(T) + F(T)) 



4 


8 


-7.524782 


-7.522773 


-7.521551 


-7.521064 


8 


29 


-7.529759 


-7.530653 


-7.530658 


-7.530674 


12 


72 


-7.530024 


-7.530175 


-7.530161 


-7.530150 


16 


145 


-7.530203 


-7.530204 


-7.530201 




18 


195 


-7.530218 


-7.530211 






20 


256 


-7.530208 


-7.530210 






24 


413 


-7.530220 


-7.530205 







k B T = 0.02E h ;E(T) 



4 


8 


-7.521660 


-7.516254 


-7.515090 


-7.514714 


8 


29 


-7.523669 


-7.524181 


-7.524207 


-7.524218 


12 


72 


-7.523698 


-7.523670 


-7.523664 


-7.523661 


16 


145 


-7.523697 


-7.523689 


-7.523687 




18 


195 


-7.523697 


-7.523701 






20 


256 


-7.523697 


-7.523701 






24 


413 


-7.523697 


-7.523698 







k B T = 0.02E h ; E(0) = ±(£(T) + F(T)) 



4 


8 


-7.529537 


-7.524058 


-7.522959 


-7.522615 


8 


29 


-7.530824 


-7.531214 


-7.531234 


-7.531242 


12 


72 


-7.530840 


-7.530824 


-7.530820 


-7.530818 


16 


145 


-7.530840 


-7.530833 


-7.530831 




18 


195 


-7.530840 


-7.530842 






20 


256 


-7.530840 


-7.530843 






24 


413 


-7.530840 


-7.530840 
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TABLE III. Ground state properties of lithium. Energies are in Eh, lattice constants in A, elastic constants in GPa. 





a 


E co h 


B 


C44 


Cn 


Cn — C12 


LDA ([3s2p]) 


3.37 


0.066 


16 


13 


18 


3.2 


LDA ([4s3p]) 


3.40 


0.066 










PWGGA ([3s2p]) 


3.44 


0.059 


14 


12 


17 


3.4 


PWGGA ([4s3p]) 


3.46 


0.059 


12 


12 


13 


3.5 


PWGGA ([4s3pld]) 


3.46 


0.059 


11 


12 


13 


3.3 


HF 


3.65 


0.020 


12 









Lit. (LDA ; PWGGA, FP-LM>W)I 
Lit. (LDA, FP-LAPW)D|-j 
Lit. (LDA, Gaussian basis )c ■ 
Lit. (LDA, plane wave basis) 
Lit. (LDA, plane wave basis) 
Lit. (LDA, Gaussian basis, KSG mode 
Lit. (LDA, Gaussian basis, HL model)l 
Lit. (LDA, Gaussian basis 
Lit. (LDA, HL, 
Lit. (MAPW. 
Lit. (HF)^ 
Lit. (HF)M 
Lit. (QMC 
Lit. (local Ansatz 
exp. 



4SK model)l 




0.062 



0.044 
0.062 
0.068 
0.061 
0.074 

0.010 (estimated HF limit: 
0.022 
0.058 
0.0517-, 
0.061E3 



0.019) 



15.0 ; 13.4 
15.4 
13.8 
13.5 
13.0 
14.7 

15.8 
14.8 
15.6 

14 
13 



17 



2.3 



13.00 11.60 14.50 2.40 



Acronyms are: 

Full-potential linear augmented plane wave (FP-LAPW) 
Rajaporal-Singhal-Kimball (RSK) 
Kohn-Sham-Gaspar (KSG) 
Hedin-Lundquist (HL) 
Quantum Monte-Carlo (QMC) 
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TABLE IV. Relative stability of the different phases of lithium. Energies are in mEh units per atom relative to the bcc 
phase, lattice constants in A, bulk moduli in GPa. In the 9R and hexagonal structures the lattice constant a refers to the 
nearest neighbour distance in the basal plane. 



phase 
fee, HF ([3s2p]) 
hep, HF ([3s2p]) 
fee, LDA ([3s2p]) 
fee, LDA ([4s3p]) 
hep, LDA ([3s2p]) 
hep, LDA ([4s3p]) 
fee, PWGGA ([3s2p]) 
fee, PWGGA ([4s3p]) 
fee, PWGGA ([4s3pld]) 
hex, PWGGA ([3s2p]) 
hex, PWGGA ([4s3p]) 
hex, PWGGA ([4s3pld]) 

fee, Ref. % LDA; PWGGA 
fee, Ref. || LDA 
hep, Ref.j], LDA 
9R, Ref. g LDA 

fee, Ref. g, FP-LAPW 
fee, Ref. 
hep, Ref. 
9R, Ref. 
fee, Ref. |8 
hex, Ref. 

fee, Ref. 

fee, Ref. 



LDA 
LDA 
LDA 
LDA 
LDA 
KSG model 
RSK model 
fee, Ref. [HI FP-LAPW 
hep, Ref.|IqFP-LAPW 

fee, RefTy, APW 
fee, Ref. |n[FP-LAPW 
fee, RefTM LMTO 
hep, Ref.j|, LMTO 
fee, Ref. Iia MAPW 



a 
4.59 
3.24 
4.25 
4.24 
3.00 
3.00 
4.33 
4.35 
4.36 
3.06 
3.08 
3.08 

4.23; 4.33 
4.28 ; 4.32 6 

3.03 
3.03 ; 3.06 6 

4.23 

4.34 

3.09 

3.07 

4.28 

3.02 

4.38 

4.20 



4.21 
4.24 



4.21 



Ebcc — E 
-0.08 a 
0.18 a 
-0.02 a 
0.19 a 
0.08 a 
0.24 a 
0.01 a 
0.04° 
0.04 a 
0.07° 
0.06 a 
0.04° 

0.15; 0.14 
0.073 ; 0.057 6 

0.062 
0.065 ; 0.050 6 

0.08 
0.09 ; 0.18 ft 

-0.01 

0.02 

0.1 

0.33 

0.25 

0.45 

0.12 

0.16 

1.4 

0.24 

0.12 

0.15 

-0.13 



B 
12 
12 
16 

16 

14 
12 
11 
14 
12 
11 

15.5; 13.3 



15.2 
13.4 
13.3 
13.3 
13.8 
13.7 
18.7 
16.8 



12.0 



a As explained in the text, the energy difference is so small that it can only be viewed as a tendency towards closed- 
packed systems. A statement about the preferred phase is not possible. 
6 Including zero point motion. 



TABLE V. Surface energies of (100), (110), and (111) surface, in units of^ 



^ 



=1.5567 



surface 

(100) 
(110) 

(111) 



LDA (3.37 A) 
[3s2p] 
0.41 
0.42 
0.49 



[3s2p] 
0.37 
0.37 
0.44 



PWGGA (3.44 
[4s3p] 
0.30 
0.32 
0.34 



[4s3pld] 
0.30 
0.32 
0.36 



Ref. ^ (LDA, 3 layers, at 3.41 A) 



0.33 
0.35 
0.40 



1G 



